Dealing with numerous reviews and widespread inquiries, it has been concluded that much more information and interpretive parameters are accessible regarding the subsurface structures when using a particular frequency range in the spectral induced polarization (SIP) measurements. Therefore, the interpretation uncertainty would diminish which causes studies with more valid and authentic outcomes. This could be achieved by using a comprehensive and general model which is appropriate for representing electrical features variation in terms of frequency, known as the Cole-Cole model. By using the SIP method and applying a defined broad of frequencies, it would be conceivable to describe items such as medium properties, spectral behavior of the studied area, and the intensity of each single parameter. The widespread use of the SIP method requires accurate and fast modeling and inversion algorithms. An integral part of every geo-electrical data inversion is an accurate and efficient forward modeling resulting in numerical simulation of responses for a given physical property model. In other words, like every other geophysical method, a reliable spectral-induced polarization inversion is highly dependent on the accuracy of the forward problem. Forward modeling is accomplished over a 2D earth structure to generate complex resistivity data by simulating current flow into the earth's surface and solving the Poisson equation containing complex values. In this contribution, a finite difference algorithm is applied to solve the complex partial differential equations (PDEs) restricted by a mixed boundary condition. A spatial Fourier transform of the PDEs, with respect to a defined range of wavenumbers, is carried out along the strike direction to elucidate 3D source characteristics. Eventually, it is necessary to conduct an inverse Fourier transform to obtain potential solutions in the spatial domain. To verify the accuracy of the proposed numerical algorithm, some synthetic models are simulated and the forward responses, including resistance and phase values with respect to a specific frequency spectrum, are calculated. Furthermore, a comparison between our numerical results and those of Geotomo geo-electrical software is provided.